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ABSTRACT 

We present a determination of the mass of the supermassive black hole (BH) and the nu- 
clear stellar orbital distribution of the elliptical galaxy Centaurus A (NGC 5128) using high- 
resolution integral-field observations of the stellar kinematics. The observations were obtained 
with SINFONI at the ESO Very Large Telescope in the near-infrared (X-band), using adap- 
tive optics to correct for the blurring effect of the earth atmosphere. The data have a spatial 
resolution of (X'17 FWHM and high S/N> 80 per spectral pixel so that the shape of the 
stellar line-of-sight velocity-distribution can be reliably extracted. We detect clear low-level 
stellar rotation, which is counter-rotating with respect to the gas. We fit axisymmetric three- 
integral dynamical models to the data to determine the best fitting values for the BH mass 
M BH = (5.5 ± 3.0) x 1O 7 M (3cr errors) and (M/L) K = (0.65 ± 0.15) in solar units. These 
values are in excellent agreement with previous determinations from the gas kinematics, and 
in particular with our own published values, extracted from the same data. This provides one 
of the cleanest gas versus stars comparisons of Mbh determination, due to the use of integral- 
field data for both dynamical tracers and due to a very well resolved BH sphere of influence 
i?BH ~ 0'.'70. We derive an accurate profile of the orbital anisotropy and we carefully test its 
reliability using spherical Jeans models with radially varying anisotropy. We find an increase in 
the tangential anisotropy close to the BH, but the spatial extent of this effect seems restricted 
to the size of i?BH instead of that Rb w 3'.'9 of the core in the surface brightness profile, 
contrary to detailed predictions of current simulations of the binary BH scouring mechanism. 
More realistic simulations would be required to draw conclusions from this observation. 

Key words: black hole physics - galaxies: individual (NGC 5128) - galaxies: elliptical and 
lenticular, cD - galaxies: kinematics and dynamics - instrumentation: adaptive optics 



1 INTRODUCTION 

The existence of supermassive black holes (BHs) in norma l galaxy 
nuclei was predicted forty years ago bv lLvnden-Belll d 1969T) . but un- 
til fifteen years ago it was still considered an interesting possibility 
which had to be demonstrated. Nowadays BHs are regarded as a 
crucial ingredient for our understanding of how galaxies form. Key 
to this paradigm shift was the launch in 1990 of the Hubble Space 
Telescope (HST). It all started with the realisation that the mass 
of the BH is correlated to other global characteristics of the host 
galaxy as a whole. Initially a correlation Mbh — L was found be- 
tween the m ass of the BH and the luminosity of the host-galaxy stel - 
lar spheroid (Kormendy & Richstone 1995; Maaorrian et al. 1998). 
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In 1997 the installation of the STIS long-slit spectrograph on HST 
allowed the spatially-resolved kinematical observations to probe in- 
side the radius of the subarcsecond BH sphere of influence Rbh = 
GMbh/o" 2 in nearby galaxies (a being the velocity dispersion of 
the stars in the galaxy). The increased accuracy in the Mbh de- 
termination s contributed to the discovery of the tighter Mb h — cr 
correlation iFerr arese & Merrill 2000; Geb hardt et al.l l2000). 

As a is a good predictor of galaxy properties, it is perhaps 
not surprising that similar correlations w ere found between M bh 
and respectively the galaxy concen tration dGraham et alfeOOll). the 



dark- halo mass dFerraresel |2002|; feaes et alj 120031; IPizzella et al 



2003 



2005), the bulge mass (McLure & Dunlop 2002; Marconi & Hunt 



Hating & Rix 2004) and the stars' gravitational binding en- 



ergy ( Aller & Richstone 2007). The existence of these correlations 
is broadly consistent with a scenario in which the BH regulates 
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the galaxy formation, during the hierarchical galaxy merging, by 
shutting off the conversion of gas into stars via a feedback mecha- 
nism due to its powerful outflows |Silk & Reesll9 98; Grana toet alj 
|2004 iDi Matteo etaDl2005l : iBower et alj2006l) . 

Even though the current scenario explains a number of ob- 
served facts, our interpretation of the role of BHs in galaxy for- 
mation is far from secure. One of the problems lies in the fact 
that the models have few observables to compare with, mainly the 
galaxy mass (or a) and Mbh. Moreover, even after a decade of 
HST spectroscopy and models, only about 4 secure BH determi- 
nations exist (compilations are given e .g. in lTremaine et al.ll2002l ; 
iFerrarese & Fordl2005l : lGraharril2008l) . With few exceptions the 
BH determinations have been performed with a single dynamical 
tracer (either gas or stars), so that no independent test of the two 
measurements methods could be made. Very little is known about 
the orbital distribution near the BHs, which i s expected to contain 
key information on the BH accretion process dOuinlan & Hernquistl 
ll997l ; lMilosavlievic & MerritfeOOll) . 

The advent of integral-field spectroscopy on all the 8-10-m 
class telescopes, combined with adaptive optics (AO) to reduce the 
blurring effect of the Earth's atmosphere is opening a new epoch 
for BH studies. The large mirrors allow for shorter exposure times 
and higher signal-to-noise ratios (S/N) of the observations, than 
what was possible with the 2.4-m mirror and the STIS long-slit 
spectrograph of HST. Observations at near-infrared wavelengths al- 
low dust absorption effects to be virtually eliminated. The integral- 
field observations, due to the tight constraint on the orbital dis- 
tribution, dramatically improve the accurac y of Mbh determina - 
tion, for a given spatial resolution and S/N dVerolme et al.|[20o3) . 
Integral-field data are also needed for a u nique recovery of the or- 
bital distribution from the observation s ( Cappellari & McDermidl 
120051 : IKrainovic et"al]l2005l ; Ivan de Ven et alj|2008l) . This can be 
understood from dimensional arguments, considering that most or- 
bits in a stationary potential conserve three isolating integrals of 
motion. This three-dimensional orbital distribution cannot be recov- 
ered without the knowledge of at least another three-dimensional 
observed quantity. Motivated by these arguments, A/bh determi- 
nations from AO-assisted integral-field obser vations in the near- 
infrared are starting to a ppear in the literature jDavies et ai]|2006l ; 
iNowak et alj2007Ll2008l) . 

The elliptical galaxy NGC 5128 (Centaurus A) is a prime can- 
didate for AO-assisted integral-field observations. At least nine in- 
dependent distance D determinations for Cen A, based on differ- 
ent methods, are available in the lite ratur f] (e.g. lTonrv et alJlioOll ; 
Reikuba 2004; Ferrare se et alj2007h . The median value is D = 3.8 
Mpc, with extreme ranges of 3.4 Mpc and 4.4 Mpc respectively. 



have been p e rformed, with ever increasing accuracy ((Marconi et alj 
200 ll 120061: iHaring-Neumaver et al.ll2006l : IKrainovic et al.ll2007k 



Here we adopt a valu4J D = 3.5 Mpc to be consisten t with all the 
earlie r papers on Mbh det e rmina ti on on this galaxy jMarconi et alj 
2001L 120061; ISilgeetal] I2005L iHlring-Neumaver et alj 12006k 
Krainovic et al.l2007l : lNeumaver e t al. 2007). At this close distance 



Cen A is the nearest elliptical galaxy and one arcsec corresponds to 
17 pc. 

Cen A is among the only ~ 10 galaxies on the whole sky with 
an observed Rbh ~ 1" (see lKormendy||2004l . for a partial list). 
Moreover Cen A possesses a nuclear gaseous disk in regular rota- 
tion from which a number of independent determinations of Mbh 



1 See the NED-ID compilation by Barry F. Madore and Ian P. Steer at 
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1 The choice of the distance D does not affect our results but sets the scale 
of our models in physical units. Specifically, lengths and masses scale as D, 
while mass-to-light ratios scale as D~ . 



Neumaver e t al. 2007). This allows for an accurate comparison be- 
tween the BH mass derived with gaseous or stellar kinematics. 
The if -band centr al surface brightnes s of Cen A is quite bright at 
fj, K « 12.2 mag djarrett et alj|2003l) . which allows a high S/N 
in the stellar spectra to be achieved in reasonable exposure times. A 
bright star, sufficiently close to the nucleus, can be used as reference 
for the AO correction. All these facts make Cen A a unique obser- 
vational benchmark for Mbh determinations in the near-infrared. 

In Section 2 we describe our data and the extraction of the stel- 
lar kinematics. In Section 3 we present the stellar dynamical mod- 
els. We discuss our results in Section 4 and we finally summarize 
them in Section 5. 



2 OBSERVATIONS AND DATA ANALYSIS 
2.1 Spectroscopic data 

For the dynamical modeling of the stellar kinematics we used 
integral-field spectroscopy obtained with SINFONI on the UT4 
(Yepun) of the Very Large Telescope of the European Southern 
Observatory on the Cerro Paranal. SINFONI con sists of the cryo- 
genic n ear-IR integral field spectrometer SPIFFI feisenhauer et alj 
2003a. bj) coupled t o the visible curvature AO system MACAO 
dBonnet et alj|2003l) . We observed the nucleus of Cen A with two 
different spatial scales: , .'250 x 07l25 (250mas scale) with a Field 
of View (FoV) of 8"x8"and 0'.'lOxO'.'05 (lOOmas scale) with a 
FoV of 3'.'2x3'.'2. The lOOmas observations a nd data reduction 
were already described in lNeumaver et al. 1 ^200% . The observations 
were performed in excellent seeing conditions of 0C'5 FWHM, as 
measured by the seeing monitor in the V-band. With the 250mas 
spatial scale the observations were taken in natural seeing, while 
with the lOOmas scale the SINFONI AO module was locked onto 
an R~14 mag star about 36"southwest of the nucleus. Although the 
reference star is relatively distant from the galaxy nucleus, the good 
seeing allowed us to achieve a nearly diffraction-limited correction 
in the if -band. Both data sets were obtained in the Jf-band, which 
covers the wavelength range 1.93-2.47/^m. The spectral resolution 
of the observations was R ~ 4800 in both scales, and corresponds 
to an instrumental dispersion of <7i ns ~ 27 km s _1 . 

The 250mas observations followed an Object-Sky-Object se- 
quence with equal exposure time of 300 s each, for a total on-source 
exposure time of 4 x 300 = 1, 200 s. The sky exposure were taken 
at 200" from the nucleus to make sure the spectra were not con- 
taminated by the large galaxy. The lOOmas observations followed 
a similar sequence with individual exposure times of 900 s each 
and a total on-source exposure time of 15 x 900 = 13, 500 s. The 
different exposures were dithered with shifts of 0>.'2 to allow for 
the removal of detector defects and cosmic rays. The data were re- 
duced using the SINFONI d ata reduction pipeline provided by ESO 
as described in section 2 of Neumaver et al. (2007). 

The spatial point-spread-function (PS F) of the lOOmas obser- 
vation was determined in section 2. 1 of Neumaver et al. I d2007l) 
from a fit to the non-thermal nucleu s, which is unresolved down to 
a / /06 (FWHM) sp atial resolution dHaring-Neumaver et alj|2006l : 
iMarconi et alj|200q) . The PSF can be approximated by two Gaus- 
sians with FWHM of 07 12 and 0730 respectively, with the smallest 
Gaussian containing 17% of the total flux. A two-dimensional fit 
with a single Gaussian to the nucleus in the reconstructed image 
from the data cube gives a FWHM of 0'.' 17. 
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For the stellar dynamical modeling the integral-field high- 
spatial resolution SINFONI observations are essential to tightly 
constrain the BH mass and stellar orbital distribution. However 
they are not sufficient as they sample only a small fr action of the 
half-l i ght radius of Cen A (Re « 83" in A-band; Ijarrett et all 
l2003h .l Shapir o et alj (2006) showed that one has to sample with 
the kinematics a significant fraction of R c to accurately constrain 
BH masses. For this reason in this paper we also use the A-band 
kinematics obtained with the Gemini Near Infr ared Spectrograph 
(GNIRS) at Gemini South bv lSilge et all d2005h . The long-slit ob- 
servations extend to a distance of R « 40" from the galaxy nucleus 
and were obtained at two position angles: (i) centered on the nu- 
cleus, along the galaxy major axis at large radii (PA p hot ~ 35° 
from North through East), and (ii) along the galaxy minor axis 
(PAphot « 125°), but offset from the nucleus by 0785. In this 
paper we increased the GNIRS a by 5% to match our SINFONI 
kinematics. Differences at this level between the two datasets are 
almost unavoidable and are likely due to low-level systematic cal- 
ibration errors, or differences in the adopted stellar templates. For 
this reason we cannot state which of the two datasets has the proper 
absolute a calibration. As none of the results in this paper depends 
on whether we apply the small shift to the GNIRS or SINFONI 
data, we arbitrarily decided to increase the GNIRS data to match 
our SINFONI ones, which we adopt as reference. A 5% uncertainty 
in a approximately translates into a 10% uncertainty in the M/L. 

2.2 Merging and binning 

The individual SINFONI exposures were registered using recon- 
structed images and merged into a single data cube. In the process 
they were resampled to a spatial scale of 0'.' 125x0'.' 125 (250mas 
scale) and Of.' 05x0'.' 05 (lOOmas scale). 

The stellar kinematics requires a high S/N for an unbiased 
extraction from the spectra. For typical observational setups, val- 
ues on the order of S/N > 50 per spectral pixel, are generally re- 
quired to detect deviations of the line -of- slight velocity distribution 
from a simple Gaussian shape ( e.g. Ivan der Marel & Franxlll993l ; 
Bend er et al ] |l994IStatleJll995l) . For this reason long-slit data are 
invariably adaptively bin ned along the spatial direction, before the 
kinematic extraction (e.g. lvan derMarelll994T) . The goal of binning 
is to maintain the maximum spatial resolution, given a constraint on 
the minimum S/N. Four our integral-field data we use the Voronoi 
binning technique|f] which constitutes an optimal so lution to this 
problem in two-dimension dCappellari & Co pin 2003). 

Given that the observations are dominated by the photon noise, 
we spatially binned the data by requiring an equal number of counts 
in each spatial bin. In this way we did not have to rely on the, often 
inaccurate, noise propagation by the reduction pipeline, to estimate 
the S/N. This constant flux per bin is achieved in the Voronoi bin- 
ning algorithm by setting the input noise N = ■s/S, where S is 
the total flux (in arbitrary units) in each unbinned spectrum. For 
the lOOmas scale the target S/N for the binning was fixed by the 
requirement for the data to be unbinned (pixel size O 7 .' 05x0'.' 05) 
within a radius R < 0'.'2 and start being binned at larger radii, where 
the surface brightness decreases. This led to 622 Voronoi bins, out 
of the original 64 2 = 4096 pixels. Similarly for the 250mas scale 
we required the data to be unbinned (pixel size 0'.'l25xO'.'l25) 
within i?<0'.'5, leading to 554 bins. The actual minimum S/N 
achieved by the adopted binning was determined afterwards, from 

3 Available from http://www-astro.physics.ox.ac.uk/~mxc/idl/ 
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Figure 1. Optimal template for Cen A. The total galaxy spectrum ob- 
served with SINFONI within the 8"x8" 250mas field (black solid line) 
is compared to the broadened optimal template determined with pPXF (red 
solid line). The fit residuals, vertically shifted by an arbitrary amount, are 
shown as blue dots. The gaps in the residuals correspond to some non well 
corrected telluric features and two H2 gas emission lines at 2.23/^m and 
2.25/^m, which were excluded from the fit. 

the rms residuals of a spectral template fit (Section [23}, to be 
(S/N)mi n ~ 80 for both spatial scales. 

2.3 Extraction of the stellar kinematics 

The near-infrared A-band spectral region is dominated by the 
strong stellar absorption feature of the 2.30/im (2-0) 12 CO band 
head. At this wavelength the galaxy spectrum is dominated by 
the light from cool and evolved giant stars. The strength of 
the CO absorption is very sensitive to the star surface gravity 
iKleinmann & Hallll98 6: Wallace & Hinkl d 19971) and this implies 
that, contrary to the optical region, in the A-band it is essen- 
tial to include giant stars in the construction of an optimal stel- 
lar template for an unb iased extraction of the stellar kinematics 
JSilge & Gebhardtll2003h . 

The extraction of the stellar kinemati cs was performed using 
the pe nalized pixel-fitting method 3 (pPXF; ICappellari & Emselleml 
2004), which fits the logarithmically-rebinned spectra with a 
template convolved with a line-of-sight v elocity-distribution de- 
scribed by a Gauss -Hermite expansion dvan der Marel & Franxl 
1 19931 ; iGerharll 19931) . The method allows the template to be care- 
fully optimized during the kinematics fit, it permits emission lines 
to be easily masked, and it includes a penalty criterion to deal with 
low S/N or insufficient resolution in the spectra. As library of stel- 
lar templates we used a set of 1 1 dwarfs and giants stars (luminosity 
class II-V) of late spectral types (K-M), observed with the same in- 
strumental setup as for the Cen A observations. 

The optimal positive linear combination of the 11 templates 
was determined only once, using a spectrum with very high S/N sa 
240, obtained by co-adding all the spectra in the 250mas SINFONI 
observations, while excluding the unresolved non-stellar nucleus. 
As expected the flux in the resulting optimal template fit is domi- 
nated by a giant M5III star (52% of the flux). This star provides by 
itself an excellent fit to the central spectrum of Cen A. However the 
contribution from a dwarf MOV star (30% of the flux) is required 
to properly reproduce the depth of the Na feature at 2.21/im (cf. 
iLvubenova e t al. 2008). The template provides a very accurate de- 
scription of the spectrum from 2.2/itn-2.4^m, which contains all 
the significant stellar absorption features in the A'-band (Fig.QJ. 

For maximum consistency the same combination of stellar 
templates was used to extract the kinematics at all spatial positions 
and for both the 250mas and lOOmas observations. Additive poly- 
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nomials of degree four were used to allow for possible variations 
in the stellar line-strength at different radii, to account for imper- 
fections in the sky subtraction or spectral calibration, and to model 
the contribution of the non-thermal nucleus (Section [2.4t . Consis- 
tent results were obtained using polynomials of degree 1-8. From 
numerous tests on other galaxies we found that, even for the most 
extreme gradients in the stellar population, the additive polynomi- 
als are sufficient to approximate possible low-frequency variations 
in the line-strength, once the global galaxy template is accurately 
determined from a high-S/iV spectrum, to account for the high- 
frequency content in the spectrum. We prefer this faster and more 
robust approach than trying to fit the template mix in every bin, 
from spectra of lower S/N. Generally both approaches give nearly 
indistinguishable results, and we tested that this is true also in the 
present case. 

We fitted with pPXF the mean velocity V, the velocity dis- 
persion a and the first two Gauss-Hermite p arameters hz and hi 
( Ivan der Marel & Franx 1993; Gerhard 1993), using the same spec- 
tral range and masked regions as in the real data. We adopted a 
penal ty A = 0.6 (defined in equation 12 of ICappellari & Emselleml 
2004). The LOSVD is always very well sampled by our observa- 
tions and the S/N is high, so any A < 0.6 gives essentially the same 
results and the measurements are virtually unbiased (Fig. |2). The 
la measurement err ors were determined as the biweight dispersion 
( Hoagl in et al. Il983l) of 100 Monte Carlo realizations. For the deter- 
mination of realistic errors we used non-penalized fits with A = 0. 
At our minimal S/N w 80 and for a ~ 150 km s _1 the typical 
random errors in V, a, hz and hi, are 4kms-\ 5 kms -1 , 0.02 
and 0.03 respectively. The kinematics extraction was performed 
while fitting precisely the same wavelength region, from 2.25/im- 
2.37/xm, in both the 250mas and lOOmas SINFONI observations 
(Fig. [5). This spectral region includes the (2-0) 1 CO band head 
and three other prominent band heads of the 12 CO series. Emission 
from the highly ionized species of [CaVIII] at 2.32p m and an H2 
line a t 2.35/nm were excluded from the fits (see also iRiffel et ail 
120081) . 

The extracted SINFONI kinematicsQiri the two scales is shown 
in Fig.[3]and Fig. [4] For the first time our observations detect a low- 
level clear sense o f stellar rotation inCen A. We used the method 
of Appendix C of IKrainovic et al. (2006) to determine the global 
kinematical major axis PA^n = 167° ± 8° and the systemic ve- 
locity of Vsyst = 531 ± 5 km s _1 . A maximum velocity of just 
AV w 20 km s _1 is reached along PAkin at about R « 4" from 
the nucleus, at the edge of our 250mas SINFONI FoV. The a field 
shows a gradual increase towards the center, before a sudden drop 
in the very center, where the non-thermal continuum prevents a 
proper kinematics extraction. The hz field is symmetric about the 
center, indicating low template-mismatch, and anti-correlated with 
V, as generally observ ed in early-type galaxies (Ben der et al J 19941 ; 
IKrainovic etatl koOS). and the hi field is generally close to zero 
over the whole field. 

The nuclear stellar rotation is counter-rotating (by about 180°) 
with respect to t he reg ular H2 nuclear gas rotation presented in 
iNeumaver et alj j2007l their fig. 6). The stars also rotate much 
slower than the gas, which reaches a maximum velocity AV ~ 130 
km s at R « / /5. This indicates that the recent gas acquisition 
was not able to produce a significant fraction of stars near the nu- 
cleus. This is consistent with the lack of evidence for any change in 
the nuclear stellar population of Cen A. This is different from a sim- 
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Figure 2. Testing of penalization in pPXF. We simulated spectra with the 
S/N = 80 of our data and an LOSVD with hz = 0.1, hi = 0.1 and 
a in the range between 30 and 240 km s _1 . We extracted the kinematics 
with pPXF and a penalty A = 0.6. The lines in the top two panels show the 
50th (median, solid line), 16th and 84th percentiles (Icr errors, dotted lines) 
of the differences between the measured values and the input values of the 
mean velocity V; n and the velocity dispersion <7; n . The bottom panels show 
the same lines for the recovered values of hz and hi, compared to the input 
values (dashed line). In the range of our observations a > 130 km s — 1 the 
recovered values are virtually unbiased. 
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Figure 3. Voronoi binned stellar kinematics from the SINFONI 250mas ob- 
servations. The four panels show the mean velocity V, the velocity disper- 
sion a and the two Gauss-Hermite moments hz and hi. Overlaid are the 
contours of the surface brightness derived from the reconstructed image, in 
1 mag intervals. North is up and East is to the left. 



ilar case of gas versus stars cou nter-rotation observed in the nucleus 
of the spiral galaxy NGC 5953 dFalcon-Barroso et alj2006t) . In that 
case the nuclear gas rotation is associated to evidence of young stars 
which co-rotate with the gas and counter-rotate with respect to the 
outer part of the galaxy. 
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Figure 4. Same as in Fig. [5] for the SINFONI lOOmas observations. The 
color levels are the same in this figure and in Fig. [5] 



2.4 Contribution from the non-thermal nucleus 

In the case of Cen A adopting a fixed template becomes impor- 
tant very close to the center, due to the presence of the non-thermal 
nucleus dMarconietafl kOQO). At radii where the nucleus domi- 
nates, the stellar absorptions are diluted, resulting in a strong de- 
crease in the observed line-strength 7. Due to the strong correla- 
tion between 7 and the stella r velocity dispersion a (section 2.2 of 
van der Marel & Franxll 19931) , an inaccurate modeling of the non- 
thermal dilution can cause large errors in the measured a. Similarly 
the rise in the non-thermal continuum could be incorrectly inter- 
preted as a variation of the stellar population, requiring a change in 
the stellar template mix and also producing an error in a. As there 
is no evidence for a sudden change in the population in the nucleus 
of Cen A the safest choice is to assume the stellar template is fixed 
and to model the non-thermal continuum via additive polynomials 
in pPXF. Additive polynomial still allow for low-order variations 
in the stellar line-strength and account for possible instrumental ef- 
fects. This approach allows the a to be reliably extracted in the 
high-resolution lOOmas observations down to 7? >0'.'2, before the 
photon noise of the nucleus eliminates all stellar information from 
the spectra. 

As an illustration of the importance of taking the non-thermal 
continuum accurately into account in the kinematic extraction, in 
Fig.[5]we show the best fit with pPXF, using a Gaussian LOSVD, to 
the spectra extracted at different radii from the lOOmas SINFONI 
data. We also performed separate fits including the Gauss-Hermite 
parameters /13 and /14. In the range / .'3-l", where we can trust 
our extraction, the LOSVD is essentially consistent with a Gaus- 
sian as /14 = —0.02 ± 0.03, in agreement with Fig. [4] The high 
S/N spectra in these plots were obtained by co-adding all the spec- 
tra within circular annuli of one pixel width (0 / .'05). The standard 
kinematics extraction approach for this paper consists of modeling 
the observed spectrum as the sum of a convolved fixed optimal tem- 
plate plus a fourth degree additive polynomial. To test the sensitiv- 
ity of the measured a to the details of the extraction, especially in 



the critical continuum-dominated nuclear region, we compared our 
standard approach with others using different options in the pPXF 
routine: (i) We use a first degree polynomial; (ii) We fit the optimal 
template at each radius using our 1 1 stars; (iii) we use as template 
a single M5III star; (iv) We include in the fit the additive contri- 
bution of a s caled version of the nuclear non-thermal spectrum (as 
described in lKelson et al.|[200oL to subtract the sky spectrum). All 
these four approaches give dispersion profiles which agree within 
the statistical errors for R > 0'.'2. 

The last approach is the one which produces the best fit and 
describes almost every details of the observations. This is the one 
we show in Fig. [5] We also give in the plots the line-strength 7, 
defined as the ratio between the mean flux within the fitted range 
in the observed Cen A spectrum, and the flux contributed by the 
best fitting stellar template alone in the same range. At J?<0 / .'2 
the non-thermal source dominates and the spectrum shows a nearly 
linear trend. At R ~ 0'.'2 it becomes possible to estimate the ve- 
locity dispersion, however the stars still contribute only 19% of the 
flux in the observed spectrum, the statistical errors are large and the 
determination is still affected by large systematic uncertainties. At 
this radius the measured a is sensitive to the details of the extrac- 
tion. Only at larger radii the kinematics can be reliably extracted 
and will be included in our models. 

Fig. [5] shows that the non-thermal continuum is still present 
in the spectra at the 25% level at R « 0'.'5. The effect of the non 
stellar continuum is still clearly visible in the varying slope of the 
spectrum. Still at all radii the stellar spectrum is well approximated 
by the fixed convolved template, once the continuum is removed 
(right column of Fig.|5j. To better quantify the extent of the dilution 
due to the central non-thermal continuum Fig. [6] shows the surface 
brightness profile I(R) measured from the reconstructed image of 
the lOOmas data cube, together with the radial profile of the the pure 
stellar contribution, estimated as E(J?) = I( R)~f(R) from the re- 
sult of the pPXF fit in the individual bins (see Ivan der Marelll 19941 
for a similar analysis on the non-thermal nucleus of M 87). The plot 
shows that the underlying galaxy profile is smooth and well approx- 
imated by a shal low power-law S(J?) oc R~ 22 , as ex pected in a 
'core' elliptical dLauer et al1ll995l ; lMarconi et al.ll200ll) . This con- 
firms that the observed radial spectral variation is consistent with 
being entirely due to the PSF effect and not to an intrinsic change 
in the stellar population. The measurement of j(R) from the spec- 
tra allows for an extremely accurate determination of the halo of 
the PSF, which would be difficult to disentangle from the underly- 
ing galaxy via photometry alone. 



3 STELLAR DYNAMICAL MODELING 
3.1 Geometry of the stellar spheroid 

Cen A appears very round on t he sky and shows little rotation 
within one R c . Wilk inson et al.l d 19861) mapped the full velocity 
field of Cen A with 71 slit positions out to 100" from the nucleus. 
They found a maximum rotation of around 40 km s " 1 roughly 
along the direction of the major axis at PAki n ~ 35°. ISilge et al.l 
(2005) find a comparable maximum rotation of about 20 km s _1 
within 4 0". Adopting a lu minosity- weighted dispersion a? « 138 
km s -1 dSilge et al]|2005l) . the amount o f rotation can be quanti- 
fied by the parameter Vmax/oo ~ 0.29 dBinnevI [19781) . From the 
empirical relation (V /a) e 0.57 (V ma x/o"o) (equation [23] of 
ICappellari et alj 12007). we estimate for Cen A (V/a) e ~ 0.17 
dBinnevI I2005IT The ellipticity of the galaxy within li? c , mea- 
sured from 2MASS if-band photometry, is very small e < 0.05 
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Figure 5. Radial variation in the spectrum of Cen A in the lOOmas SIN- 
FONI observations. Left Column: Different panels show the observed spec- 
tra (black solid line) obtained by co-adding the spectra of the spaxels con- 
tained within circular annuli of radius R and one pixel width (f/.'05). The 
best-fitting pPXF model (red solid line) consists of the stellar template 
(Fig. [TJ convolved with a Gaussian LOSVD plus a fourth degree additive 
polynomial, plus a scaled copy of the non-thermal nuclear spectrum (top 
panel). The residuals are shown at the bottom of each panel with the blue 
dots. The value 7 > 100% at R = 1" is likely due to a small difference 
in the line strength of the stellar spectrum, or in the sky subtraction, with 
respect to the adopted optimal template. Emission from the [CaVIII] line at 
2.32/xm and the H2 line at 2.35/im were excluded from the fits (grey areas). 
These lines are prominent also outside the nucleus and cannot be properly 
subtracted by our model. Right Column: The convolved optimal template 
(red solid line) is compared to the observed spectrum after subtraction of 
the nuclear spectrum and polynomial contributions (black solid line). The 
blue dots show the residuals. 
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Figure 6. The radial surface brightness profile I(R) as measured from the 
individual pixels in the SINFONI lOOmas scale (filled squares) is compared 
to the estimated profile of the stellar light S(i?) = I(R)f(R) (red open 
circles). As expected for a normal galaxy profile, the underlying stellar light 
profile is smooth and well approximated by a power-law <x R~ 0,22 

(solid line). The dashed vertical line marks the radius where the PSF has 
reached 50% of its peak value. 



dSilgeetal.ll2005l) . These measured values seem to place Cen A 
on the region of the face-on fa st rotator early-type g alaxies on the 
(V/a, e) diagram (figure 1 1 of ICa ppellari et al . 2007). 

Although the ellipticity of Cen A remains low e < 0.1 out to 
about 3R C , a well-defined photometric PA can still be determined. 
The major axis that we measure from the inertia elli psoid of the 
surfac e brightness on the 2MASS 7f-band photometry jjarrett et alj 
120031) at large radii is PA phot = 35° ± 3° (East of North). The axis 
of maximum rotation that we determine from our 250mas SINFONI 
data is quite different at PA ki „ = 167° ± 8° f Section |2"73||. The 
sense of stella r rotation that w e meas ure is in support to the de- 
termination by Wilkinson et al. ( 1986), and our accurate PA k m is 
probably consistent with their determination, given their large un- 
certainties. 

The strong kinematical misalignment of ~ 48° is an indica- 
tion that the galaxy is certainly not an axisymmetric object. This 
is consistent with the strong twist in the kinemat ical position an- 
gle of the planeta ry nebula system at large radii dHui et alii 19951 : 
IPeng et al] |2004b). However the galaxy is unlikely to be well de- 
scribed by a stationary triaxial geometry either. In fact at radii where 
Cen A starts to appear more elongated, its surface brightness is 
dominated by stellar shells, which constitute the relics of a recent 
accretion event and are likely associated with the polar geometry of 
the st r ong dust lane which crosses the g alaxy nucleus dMalin et"aT] 
Il983l ; lOuillen et alj|l99j llsraell 1 19981) . This shows that at those 
radii the galaxy has not reached an equilibrium configuration. The 
relatively low but ordered rotation near the nucleus, the presence 
of prominent stellar shells at large radii and the strong kinemati- 
cal misalign ment of Cen A are rem iniscent of the fast rotator SO 
NGC 474 in lEmsellem elail fe007h . However only much more ex- 
tended integral-field observations in the near-infrared, sampling up 
to IR C , could conclusively reveal the dynamical status of the central 
regions of Cen A. 

At different radii the galaxy is expected to reach an equi- 
librium configuration in very different time scales, starting from 
the center. We estimated the circular v elocity V c of Cen A 
from the Multi-Gaussian Expansion (MGE uEmsellem et al. I I 19941) 
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parametrization of its if -b and surface brightness tabulated in 
Haring-Neumaver et al. U20061) and the best fitting mass-to-light ra- 
tio M/L we derive in this paper. It is nearly flat over two orders of 
magnitude in radius V c ~ 200 km s _1 from 4"^100". This implies 
that at the edge of our SINFONI observations (R ~ 4") the char- 
acteristic orbital period is T ~ 2 Myr, while it is T ~ 0.2 Gyr at 
large radii where the stellar body starts to become more elongated 
and shells are still visible. The very short orbital period in the cen- 
tral regions, where the galaxy is nearly circular in projection, sug- 
gests that it has already reached an equilibrium configuration there. 
Given that the kinematics at a certain projected radius R is mostly 
influenced by stellar or bits having that same characteristic radius R 
dKrainovic et alj|2005h . this justifies the use of stationary dynami- 
cal models to study the central dynamics of Cen A, even though the 
galaxy is still dynamically evolving at larger radii. 

For this paper we will construct two types of models, to test 
the sensitivity of the A/bh estimate to the assumed geometry and 
dynamics. We build (i) an axisymmetric orbit-based model, to re- 
produce in detail the kinematic observations, and (ii) a simple 
anisotropic spherical Jeans model to qualitatively check the Mbh 
determination. 



3.2 Axisymmetric three-integral models 

Our axisymmetric th ree-integral dynamical model is based on 
ISchwarzschildl d 19791) numerical orbit-superposition method. This 
has become the current standard for all Mbh determinations 
from the stellar k inematics, available in the literature from dif- 
ferent groups (e.g. Ivan der Marel et al.lfl998l : iGebhardt et~al]|2003l ; 
IValluri et alj|200l) . The axisymmetric implementation we use, the 
orbital 'dithering' and the se tup we employ, are described in de- 
tail in lCappellari et ID §006). The relatively limited spatial exten- 
sion of the integral-field data, and the likely lack of equilibrium at 
large radii, does not justify the use of more genera l triaxial models 
dde Lorenzi et alj|2007l;lvan den Bosch et al.ll2008b . 

Our dynamical model assumes constant M/L. This is likely 
a good approximation in the central regions (R < R a ) of Cen A, 
where dark matter is expected to contribute only a small fraction 
of the mass etajj 1995; Peng et al. 2004b) as observed via 



dynamics or gravitational len s ing in larger samples of early-type 

Cappellari et al 



galax ies i Gerhard et al. 
2006 ; iKoopmans et al.l 



__np 

2001; iRusin et alj 120031; ICati.. 
20061; iThomas et alj 120071 ; iBolton et al 



2008). The model still allows for dark matter in the form of a 



constant shift in the global M / L. Our model assumes an MGE 
parametrization for th e jf-band surface brightness of Cen A using 
the values tabulated in ;Haring-Neumaver et alJd2006l) . 

The models are fitted to our integral-field SINFONI lOOmas 
(with AO correction) and 250mas kinematics in the cen tral regions 
(R<4 "), and to the two-slits GNIRS kinematics of Isilge et ail 
d2005h at larger radii (R < 40"). For R < 1" we only used the SIN- 
FONI lOOmas data in the fit, not to spoil the information contained 
in the high-resolution AO-assisted observations with the seeing- 
limited 250mas ones. Similarly we do not fit the GNIRS kinematics 
for R < 3'.'5, where we have higher quality SINFONI data. As the 
models are bi-symmetric by construction, the SINFONI kinemat- 
ics were bi-symmetrized before the fit along PAkm = 167°. For 
all datasets we fit, and we compute the predictions, for the Gauss- 
Hermite moments up to h^-he, where h$ and hi are the measured 
values, while we set hs — h?, — ± 0.3 (la errors). This is done 
to include in the model the additional observational constraint that 
extreme values for the high-order moments are never ob served (see 
two extreme examples in fig. 4 of ICappellari et alJ 2007). 



We constructed a set of models with slightly different ge- 
ometries to test the variation in the recovered Mbh'- (i) An edge- 
on model (i = 90°); (ii) An edge-on model with the major axis 
GNIRS kinematics aligned with PAkm = 167°. In this case we 
set V — h-j, — for the GNIRS minor axis data; (iii) An edge-on 
model with 'circularized' GNIRS kinematics. Here we assumed the 
a and /14 field have circular symmetry, and a profile given by the 
average of the two slits. We fitted the a and /14 GNIRS kinematics 
along five polar sectors, covering the whole galaxy from the ma- 
jor to the minor axis, and not just the observed major/minor axes, 
(iv) Same as model [ii], but with i = 60°, which is the average in- 
clination for random orientations; (v) Same as model [ii], but with 
i = 45° . All these models are expected to produce similar results 
given that the differences between the major and minor axis GNIRS 
kinematics are almost at the level of the measurement errors, and 
given that Cen A is so close to circular within the region where 
we have kinematics, that for large ranges of inclination the intrinsic 
density is essentially spherical. We find that indeed all these models 
produce BH mass estimates well within the errors, and have simi- 
lar uncertainties. The orbital distribution of all the models is also 
nearly the same. For this reason we will not present the individual 
results. Instead for all our figures below we adopt as representative 
the model (ii). 

The two best-fitting parameters of the model, A/bh and M/L, 
were determined as follows. After an initial approximate search 
for the best fitting parameters using isotropic Jeans models (Sec- 
tion [33), we computed a set of orbit libraries, each consisting of 
444,528 individual orbits, at our first guess for M / L. For each set 
the Mbh was sampled linearly in u — sinh(MBH /A/bh, best), 
where A/bh, best is our first guess for A/bh- We then scaled 
the velocities of the orbit s to compute models at different M/L 
dvan der Marel et alJ,ll998l) . In the fit we use a mode st amount of 
regularization A = 10 (see Ivan der Marel et all l 19981 , for a defini- 
tion). The contours of the resulting x 2 , which describe the agree- 
ment between the models and the data, are shown in Fig. UJ The 
best fitting BH has mass A/bh = (5.5 ± 3.0) x 10 7 M Q , where 
the 3<r error bars are marginalized over the M/ L and correspond to 
one degree-of-freedom (A% 2 = 9). The corresponding best fitting 
/C-band M/L = (0.65 ± 0.15) M /L . 

The data-model comparison for our best fitting model are 
shown in Fig. [8] All the significant features of the data are repro- 
duced in detail by the model, for the two SINFONI datasets. Unsur- 
prisingly the GNIRS long-slit data are also reproduced in detail (not 
shown). Although we only have observations along two slits, we 
still predict the model kinematics over a full quadrant, to verify that 
the model does not contain un-physical features o r sharp disconti- 
nuitie s, which are never observed in real galaxies dEmsellem et alj 
120041) . 

The orbital anisotropy of our best fitting model is shown in 
Fig. [9] We plot the ratio a r /a t between the radial velocity dis- 
persion oy and the tangential one defined as a\ = [cA + cr|)/2. 
Here 0$ includes only random motion and not ordered rotation, so 
that an isotropic system has a r /a t = 1, and (r, 6, <j)) are the stan- 
dard spherical coordinates. We also plot for reference the radius at 
which the total luminous mass in the MGE model is equal to A/bh, 
for our best fitting M/L and A/bh values. The derived value is 
7?bh = 0'.'70, which is very close to the usual practical definition 
of the BH radius of influenc e Rbb = GMbb/ oI = 0'.'68, for 
a e = 138 kms -1 (see also lMarconi et al.l 2006). Also shown is 
the break radius R b = 3'.' 9 of a Nuker-law dLauer et al.ll 1995b fit 
to the /("-band HST photometry dMarconi et alj|2000h - Both radii 
are very well resolved by the observations, which then allow the 



© 2008 RAS, MNRAS 000,[QfT2] 



8 Cappellari et al. 




2.0 



Figure 7. Determination of the best fitting parameters. The contours of the 
X 2 , which describes the agreement between the data and the model, are 
plotted as a function of Mbh and X-band M/L. The lowest levels show 



the Ax 



X~ 



Xn 



which correspond to confidence levels of 1 , 2 and 



3<t (thick contour), for one parameter. Additional contours are separated 
by factors of two in Ax 2 . The o pen square marks the b est fitting Mbh 
derived from the gas kinematics bv lNeumaver et alj 420071) . The open circle 
is the best fitting isotropic model from Section 13.31 The small dots are the 
locations of the models that were run, while the x 2 was interpolated at other 
locations. 
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Figure 8. Data-model comparison for the best fitting three integral model. 
Top two panels: The top row shows the symmetrized and linearly interpo- 
lated lOOmas SINFONI data of Fig.g] The second row show the best fitting 
dynamical model predictions. The central bins that were excluded from the 
fit are shown with the white diamonds. Bottom two panels: Same as in the 
top two panels, for the 250mas SINFONI kinematics. For each quantity the 
color scale is the same in the two instrumental configurations. 




Figure 9. Anisotropy variation. The diamonds, connected with solid lines, 
show the anisotropy o> / at (see text for a definition) measured at different 
polar angles, from the equatorial plane to the symmetry axis, in the galaxy 
meridional plane. As the model is nearly spherical, the differences at a given 
radius, provide a rough indication of the model uncertainties. Also shown 
are the break radius (or core radius) Rf,, the BH radius of influence Rbh 
and the a of the best fitting Gaussian model of the SINFONI 1 OOmas PSF. 



nuclear orbital distribution to be robustly recovered. The velocity 
ellipsoid appears to deviate from a sphere (isotropy) by just 15% 
in the whole range where we have ki nematics. This is just above 
the errors of 5% in this determination ( Cappellari et al. 2007]). The 
main significant trend is a decrease in the a r /a t ratio inside -Rbh, 
but no sharp transition is seen near Rb- 



3.3 Anisotropic spherical Jeans models 

In the previous section we used general axisymmetric Schwarz- 
schild models to determine Mbh and to derive the anisotropy 
profile in the nucleus of Cen A. In this section we use simpler 
anisotropic spherical Jeans models to test the reliability of these 
modeling results. 

In the past 25 years the dynamical models of galaxies have 
evolved from simple sp herical models based on the Jeans equations 
teinnev & Mamonlll982T) to quite g eneral orbit- or particle- based 
numerica l techn iques with spherical dRichstone & Tremaindl 19881 : 
iRix et al.1 Il997h . axi symmetric (e.g. [van der Marel et al. 1998) 
or tri axial geometry ide Lorenzi et alj 20071 ; Ivan den Bosch et al. 
2008). Contrary to the models based on the Jeans equations, which 
can generally be solved and tested to machine precision, the results 
of the more general models are sensitive to the numerical imple- 
mentation details and more difficult to test thoroughly. Moreover 
the complexity of the general models often prevents a qualitative 
assessment of the results, which sometimes h ave to rely entirely on 
subtle detai ls of the kinemati cs (e.g. fig. 2 of iGebhardt et alj|2003l 
or fig. 7 of lNowaketafll2007l) . 

One way to gain confidence in the A/bh determinations and 
the recovered anisotropy is by comparing them with the simpler 
Jeans models ( Ivan der Marel et al]|l998l : ICretton & van den Boschl 
1999). These models are based on a completely different set of as- 
sumptions and do not suffer from possible incompleteness in the 
sampling of the orbit library or from limited kinematic coverage. 
To test our best fitting A/bh and M/L and the corresponding re- 
covered anisotropy profile of Fig. [9] we use the Jeans Anisotropic 
MGE (JAM) package 3 of Cappellari (2008). The method allows 
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one to compute the predicted second moments, projected onto the 
sky plane, for a model with a variable anisotropy profile, via a quick 
and accurate single numerical quadrature. Given that our standard 
model of Cen A is intrinsically cl ose to spherical, w e use the spher- 
ical formalism of equation (50) o f lCappellarjj2008t) . 

We construct a spherical model having the same MGE surface 
brightness as our Schwarzschild model, but we set the axial ratios 
of all the Gaussians q[ — 1. We assume the same anisotropy profile 
of our best fitting Schwarzschild model of Section \T2\ with the 
same Mbh and M/L. The method allows for a different anisotropy 
P = 1 — <r| /a 2 to be assigned to different Gaussian components 
of the MGE model. In practice we assign an anisotropy = 0.2 
(radial anisotropy) to the Gaussians with 0'.'5 < <Ji < 10" and 
an anisotropy j3 — —0.35 (tangential anisotropy) to the remaining 
Gaussians. Although /3 changes in a discontinuous fashion for the 
different MGE Gaussians, the resulting P profile is smooth. In fact 
the Gaussians overlap with each other and the P at a certain spatial 
position is a luminosity-weighted sum of the /3 of the individual 
components. 

This anisotropic model, mimicking the best-fitting three- 
integral model, is shown with the solid line in the top panel of 
Fig. [10] Although the model is not by itself a fit to the data, it 
provides a remarkably accurate description of the observed and 
radially-averaged SINFONI and GNIRS (outside R > 1") second 
velocity moments V TUla = \/V 2 + a 2 . We ignore the h,4 Gauss- 
Hermite moment in the Kms estimation, as it is essentially zero 
within the errors. The purely Gaussian LOSVD extraction of Fig. [5] 
gives consistent Kms values within the errors. This test shows that 
the best-fitting discrete Schwarzschild's numerical representation of 
Cen A is an excellent approximation for the continuum limit repre- 
sented by the Jeans model. It provides a strong confirmation for the 
recovered anisotropy profile and the need for tangential anisotropy 
in the nucleus of Cen A. We also show in Fig. [TJ5J with the two 
dashed lines, the Jeans predictions for two Mbh corresponding 
to the 3cr upper and lower confidence limits of the Schwarzschild 
model. They show that, as expected, inside R < 1" the data are very 
sensitive to a change in Mbh at the level of the quoted errors. 

The dotted line shows the model prediction, still with the same 
P profile as for the previous models, fo r Mbh = 2 x 10 s Mp. This 
is the best fitting value determined by ISilge et alj (2005) from the 
GNIRS data. This model has a much steeper Kms profile near the 
center, which seems to qualitatively reproduce the steeper rise of 
the Vrms in the GNIRS data within R < 1. Unfortunately we cannot 
compare our Jeans models with the more nuclear GNIRS kinemat- 
ics (R <0'.'5). Those measurements have rather extreme hi >0.2 
values. In this case it becomes not possible to reliably translate 
the kinematics into a true Kms value, as required by the Jeans 
equations. In fact when the Gauss-Hermite moments are large, the 
Vrms one can derive by formally integrating over the LOSVD is 
extremely sensitive to the wings of the LOSVD which are observa- 
tionally not well constrained. For this reason we limit our compar- 
ison to the remaining values, for which hi is consistent with zero 
within the lcr errors. Our comparison suggests that the ~ 4x higher 
value for Mbh derived by those authors is not due to differences in 
the modeling, but more likely to the differences in the measured nu- 
clear kinematics profile with GNIRS. A more conclusive evidence 
would require a full Schwarzschild model of the nuclear GNIRS 
data, which is outside the scope of this paper. 

It is clear however that our S INFONI nuclear ki nematics is not 
consistent with the GNIRS one of lSilge et all J2005h when R < 1". 
Those authors measure larger a values and especially strong hi val- 
ues in the nucleus, while we find hi ~ over the full field-of-view 
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Figure 10. Spherical Jeans models. Top Panel: Our SINFONI observations 
of V rms = VV 2 + a 2 for the lOOmas (red open circles) and 250mas scale 
(blue open squares), biweight averaged over circular annuli, are plotted as 
a function of the radius from the center of Cen A. This is the same kine- 
matics presented in Fig. 1 3 141 which was used in the Schwarzschild mod- 
els. We also show for comparison, as black open upward triangles, the a 
measured in Fig. fj] assuming a pu rely Gaussi an LOSVD. The correspond- 
ing GNIRS Kms observations o f ISilge et alj koOih along the major axis 
(filled downward triangle) and the minor axis (filled upward triangle), are 
also shown for R > 0'.' 5 (see text for a discussion of the discrepancy be- 
tween the GNIRS and SINFONI data inside R < 1"). These have been ar- 
bitrarily scaled up by 5% to match our data at R > 1" (see text). The solid 
line shows the prediction of a JAM model, having the best fitting M^n 
and M/L derived from the Schwarzschild model of Fig. [7] and the corre- 
sponding anisotropy profile of Fig. [9] The model was convolved with the 
PSF of the lOOmas observations. The two dashed lines have a A/bh corre- 
sponding to the upper and lower 3cr errors in the Schwarzschild model. The 
d otted line has Mb h = 2 x 10 s , corresponding to the best fitting value 
of lSilgeetal3 l2005). Bottom Panel: Same as in the top panel, but assuming 
isotropy (j9 = 0) for all the Jeans models. The M/L has been decreased to 
match the data at large radii. 



(Fig.|4j. The reasons for these differences are likely the following: 
(i) Our SINFONI data represent a significant improvement over the 
GNIR S one, both i n S/N and in spatial resolution (compare fig. 6 
o flSilge et alj2005l with the right column of our Fig.^. The GNIRS 
data forced them to extract the kinematics inside R < 1" by fitting 
only the first CO band, while we could accurately fit all three CO 
bands from the SINFONI data, down to the limit where they disap- 
pear due to the rise of the non-thermal nucleus; (ii) We used differ- 
ent approaches for the subtraction of the non-thermal continuum. 
As discussed in detail in Section [2~4l a non accurate subtraction of 
the non-thermal contribution from the spectrum, leads directly to 
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an error in a. We subtracted the non-thermal source with additive 
polynomials during the pPXF fit, so that both the kinematics and 
the continu um contribution are simultaneously optimized to fit the 
spectrum. ISilgeetal .1 (2005) instead removed the continuum based 
on a measure of CO line-strength and only subsequently fitted the 
kinematics on the continuum-free spectrum. 

For comparison, in the bottom panel of Fig. [To] we show four 
spherical Jeans models with the same A/bh as in the top panel, but 
in this case assuming isotropy {J3 = 0). The M/L was decreased 
to fit the data at large radii, where the BH contribution is minimal. 
A reasonable fit to the Kms is obtained in this case at the upper 
limit M BH = 8 x 10 7 M Q of our Schwarzschild confidence inter- 
val. This general agreement confirms previous comparisons, show- 
ing that although general three-integral models are needed for accu- 
rate Mbh determinations, much simpler isotropic models, applied 
to the same data, provide a useful sanity-check for the determi- 
nation (van der Marel et alj|l998l: ICretton & van den Bosch l ll999l : 
Ijoseph et alJl200ll ; IVerolme et alil2002h . These simple tests would 
be especially useful to give confidence in the A/bh determinations 
in critical cases, where Rbh is barely resolved and a determination 
of A/bh see ms at the limit of what is possible with curr ent tech- 
nolo gies (e.g.lDavies et alJl2006l : lNowak et al .|[2007l . 1 200 8h . 



Gebhardt (2004) compared the A/bh determinations between 



two-integral and three-integral models for 10 galaxies and found 
that, eve n using seeing-limited ground-based data, the two-integral 
Mbh of lMagorrian et al.l d 199a) are systematically too high by just 
a factor ~ 2, with respect to more recent three-integral models, 
based on high-resolution HST observations, ft would be interesting 
to perform a similar type of comparison, but using the same data for 
both anisotropic Jeans and Schwarzschild models. 

The reason for the reasonable agreement between the simpler 
and the more general models is likely due to the fact that many real 
galaxies do not 'utilize' all the freedom in the orbital distribution 
that would be required to produce more dramatic disagreements, 
and do n ot deviate very strongly from being isotrop i c in their cen- 
ters ( e.g. IVerolme et alj|2002K iGebhardt et a l. 2003; Shapi ro et alj 
2006). This is certainly true for the results presented in this paper. 



4 DISCUSSION 

Our determination, using the stellar kinematics, for the BH mass 
in the center of Cen A A/bh = (5.5 ± 3.0) x 1O 7 M agrees 
very well with the determination Mbh = (4.9 ± 1.4) x 10 7 M Q 
we performed using the gas kinematics in iNeumaver et alj d2007l) 
from the same SINFONI data. This gas-stars A/bh comparison con- 
stitutes one of the most robust and accurate ones, due to a very 
well resolved BH sphere of influence (Rbh ~ 0'.'70 compared to 
a PSF FWHIVfe 0'.' 17) and thanks to the use of high-resolution 
integral-field data for both kinematical trac ers. Comparable and 
equ ally successful com parisons were done bv lShapiro et al.l (2006) 
and Siopis et al. ( 2008), while a less good agreement was found in 
Capp eilari et alj (120021), li k ely du e to the disturbed gas kinematics. 

In INeumaver et alj 1(2007) we made a detailed com- 
parison with all the numerous previous Mbh determina- 



tions for Cen A dMareoni et all 1200 1[ 12004 Isilee et al.l [20051 : 
Haring-Neumav er et alj|2006l : iKrainovic et alj|2007l) . In summary 
the only significant disagree ment in Mbh is w ith the previous stel- 
lar Mbh determination by Sil ge et all d2005l) . In Section 13.31 we 
showed that the disagreement is likely caused by a difference in the 
data quality and in the treatment of the contribution from the cen- 
tral non-thermal continuum in the kinematic extraction. It is not due 



to differences in the de tails of the modeling methods. As shown in 
INeumaver et alj d2007t) . with our new gas and stars Mbh determi- 
nation C en A lies within t he erro rs on the Mbh — o relation as given 
by either iTremaine et"al] d2002l) or lFerrarese & Fordl d2005l) . 

The value (M/ L)k = 0.65 is well consistent with what one 
would expect from a ~ 7 Gyr luminosity-weighted age of the stellar 
population of Cen A, almost ind ependen t ly of t he assumed near- 
solar metallicity. This adopts the iKroupal d200lh normalization of 
the initial mass function (IMF) and uses the population models of 
iMarastorJ (2005). These population models include a proper treat- 
ment of the TP-AGB stars which is essential for a reliable prediction 
of the near-infrared flux and (M/L)k- The (M/L)k of the stellar 
population would decrease to 50% of the observed dynamical value 
if Cen A had a mean age of ~ 3 Gyr. In that case dark-matter would 
be required to explain the observations. However the younger age 
is at the lower extreme of the estimated age of 3-8 Gyr, for the 
mai n body of Cen A, f rom the analysis of its globular cl uster sys- 
tem l Peng et al.l2004al) . Adopting the normalization of the[Salpeter 
JT955F IMF, the mean age of Cen A would have to be lower than 
~ 3 Gyr for the population (M/L)k not to over-predict the dy- 
namical one. However this IMF normalization has been shown to 
be inconsistent wi th the observed dynam ica l M/L of early-type 
galaxies as a class dCappellari et ai]|2006l : see|de Jong & Bel1l2007l 
for a review). 

Cen A possesses a clear core in the luminosity profile, 
with a break from a shallow cen tral slope to a steep er outer 
one at a radius R b ^ 3'.' 9 dMarconi et alj |2000|) , which 
is extremely well resolved by our kine matics. Every galaxy 
spheroid seems to host a supermassive BH dMagorrian et alJI 19981 ; 



ati (IMago 

iFerrarese & Merrittll200b1 ; IGebhardt et alj|2000l) . and galaxies are 
expected to form by mergers. These cores have been interpreted 
as due to the scouring of the galaxy profile due to the ejection 
of stars in radial orbits p assing close to the BH binary which 
forms shortly after mergers 1 Faber et aljl997l : l0uin lan & Hernquist 



ll997l : lMilosavlievic & Merrittll200ll) . The models predicts that the 
anisotropy should show a bias towards tangential orbits inside Rt- 

To test this prediction, in this paper we present the recovered 
stellar anisotropy profile o r /a t for Cen A (Fig. [9). The profile is 
to first order flat around the value o r /o t ~ 1. The maximum de- 
viations of the velocity ellipsoid from a sphere (isotropy) are on 
the order of 15%. The profile shows no sudden bias towards tan- 
gential orbits inside Rb, while a tangential bias appears only in- 
side Rbh- This profile seems typical of the few nuclear anisotropy 
profiles that have been published so far using integral-field observa- 
tions or spherical models (the only cases for which the orbital distri- 
bution can be robustly recovered ): M32 dVerolme et al.l2002 ), M8 7 
dCappellari & McDermidll2005r). NGC 3379 dShapiro et al.N2006h . 
NGC 1399 dHoughton et alj|2006l: IGebhardt et alj|2007l) . A similar 
anisotr opy trend was observ ed with high-resolution HST long-slit 
data bv lGebhardt et al. although in that case the profiles ap- 

pear more noisy. 

Although the predicted tangential bias is observed, its size 
does not seem related to that of the core. So these anisotropy de- 
terminations would not seem to support the scenario in which the 
cores have been produced by the core-scouring mechanism. How- 
ever a big caveat is that the current simulations are not entirely re- 
alistic: they only consider binary mergers of cusp galaxies which 
are initially isotropic. Now that high-quality anisotropy measure- 
ments are starting to appear, it would be valuable if the numerical 
simulations could revisit their predictions in a more cosmological 
motivated setting, considering sequences of multiple mergers. For 
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the moment we must refrain from making firm statements on what 
our anisotropy profile implies for BH formation. 



5 SUMMARY 

We present a determination of the mass of the supermassive BH in 
the nucleus of the nearby elliptical galaxy Cen A, using high spa- 
tial resolution integral-field observations of the stellar kinematics, 
combined with large-scale long-slit kinematics. Our high-resolution 
observations were obtained in the near infrared i^-band with SIN- 
FONI on the VLT, using adaptive optics to correct for the blur- 
ring effect of the atmosphere. The PSF of our data has a close 
to diffraction-limited resolution of 0717 FWHM and very high 
S/N > 80 per spectral pixel. We discuss the extraction of the stellar 
kinematics, with particular emphasis on the treatment of the non- 
thermal continuum coming from the nucleus of Cen A, which con- 
taminates the central stellar spectra. 

We use a standard three-integral axisymmetric numerical 
orbit-superposition model to determine the mass of the supermas- 
sive BH, which turns out to be M BH = (5.5 ± 3.0) x 1O 7 M . This 
value is in good agreement with our previous accurate determina- 
tion Mbh = (4.9 ± 1.4) x 10 7 M Q , obtained fr om the gas kine- 
matic s extracted from the same SINFONI data iNeuma ver et all 
2007). It is also i n agreement w i th the recent analysis of the HST 
gas kinematics bv lMarconi et al. Illooi), and with other similar, but 
with lower spatial resolution, Mbh determinations. This mass is 
consistent with the pr edic tion of the Mbh — a re lation of either 
Tremai ne et al. | d2002h or iFerrarese & Fordl d2005h . We carefully 
test the Mbh and anisotropy recovery from our three-integral stel- 
lar modeling using simple anisotropic spherical MGE Jeans mod- 
els. We find very good consistency between the two different ap- 
proaches. 

We study the nuclear orbital distribution of Cen A. In agree- 
ment with a few previous studies of other galaxies, we find a tan- 
gential bias in the anisotropy near the BH. However the size of the 
tangential region do not seem to be associated with that of the core 
radius Rt of the surface brightness profile, but has the size of the 
much smaller BH sphere of influence -Rbh- This does not add sup- 
port to the scenario in which galaxy cores are scoured by binary 
BHs. More realistic numerical simulations would be required for 
a comparison with the new accurate data and to draw conclusions 
from this observation. 
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